What we want:
In order to understand the $\chi^2$ test used a toy a model. Remembering that the estimator is:
After that we used the $\Delta \chi^2$ to find the confidence levels of our parameters.
In [1]:
import sys, platform, os
import matplotlib
matplotlib.use("nbagg")
from matplotlib import pyplot as plt
import numpy as np
from matplotlib import rc
rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
## for Palatino and other serif fonts use:
#rc('font',**{'family':'serif','serif':['Palatino']})
rc('text', usetex=True)
#LaTeX
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
#uncomment this if you are running remotely and want to keep in synch with repo changes
#if platform.system()!='Windows':
# !cd $HOME/git/camb; git pull github master; git log -1
print('Using CAMB installed at '+ os.path.realpath(os.path.join(os.getcwd(),'..')))
sys.path.insert(0,os.path.realpath(os.path.join(os.getcwd(),'..')))
import camb
from camb import model, initialpower
In [2]:
import scipy.integrate as integrate
import scipy.special as special
from scipy.interpolate import interp1d
# Anderson et al. 2013
tab_cov_xi0 =np.loadtxt('tab_cov_xi0_DR10.txt')
tab_bins = np.loadtxt('tab_bins_DR10.txt')
cov_mat = np.matrix(tab_cov_xi0)
#print cov_mat.I
cov_inv = cov_mat.I
nbb = len(tab_bins)
tab_err_xi0 = np.zeros(nbb)
In [ ]:
plt.matshow(cov_mat)
plt.show()
plt.matshow(cov_inv)
plt.show()
In [3]:
# Link between the distribution of total matter and our tracers
#(Galaxies). Want to constrain cold dark matter distribution.
nb_bias = 100
tab_bias=np.linspace(1., 3., nb_bias)
#
arrob = np.linspace(0.009,0.050,30)
#
tab_chi2 = np.zeros((len(arrob),nb_bias))
tab_chi2_cov = np.zeros((len(arrob),nb_bias))
In [4]:
r, xi0, xi2 = np.loadtxt('Anderson_2013_CMASSDR10_corfun_x0x2_postrecon.txt')
#Covariance matrix start at 46 Mpc/h
rcut = 46.
ind_cut = int(np.argmin(abs(r-rcut)))
nbb = len(r[ind_cut :])
tab_bins = r[ind_cut :]
tab_xi0 = xi0[ind_cut :]
tab_xi2 = xi2[ind_cut :]
In [5]:
for indx in range(len(arrob)):
pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=arrob[indx], omch2=0.122)
pars.set_dark_energy() #re-set defaults
pars.InitPower.set_params(ns=0.965)
#Not non-linear corrections couples to smaller scales than you want
pars.set_matter_power(redshifts=[ 0.55], kmax=10.0)
#Linear spectra
#pars.NonLinear = model.NonLinear_none
#Non-Linear spectra
pars.NonLinear = model.NonLinear_none
results = camb.get_results(pars)
kh, z, pk = results.get_matter_power_spectrum(minkh=1e-4, maxkh=10, npoints = 5000)
s8 = np.array(results.get_sigma8())
xi = np.zeros(nbb)
for i in range(0,nbb):
tab_err_xi0[i] = tab_cov_xi0[i,i]
kstart = 0
kcut = 10
kcut_ind = np.argmin( abs(kh - kcut))
factor = np.power(kh[kstart:kcut_ind],2) * pk[0,kstart:kcut_ind] / (kh[kstart:kcut_ind]*2*np.pi**2)
for i in range(0, nbb):
IntegrandXi0 = factor * np.sin(kh[kstart:kcut_ind]*tab_bins[i]) / tab_bins[i]
# 2pt correlation function
xi[i] = np.trapz(IntegrandXi0,kh[kstart:kcut_ind])
vals_err = np.power(tab_bins,2)*np.sqrt(tab_err_xi0)
vals_measure = np.power(tab_bins,2)*tab_xi0
for i in range(0, nb_bias):
vals_mod = np.power(tab_bins,2)*xi*tab_bias[i]**2
vals_mod_cov = xi*tab_bias[i]**2
tab_chi2[indx,i] = np.sum( np.power((vals_mod - vals_measure)/vals_err,2) )
tab_chi2_cov[indx,i] = np.dot( vals_mod_cov-tab_xi0 , np.dot( cov_inv ,
(vals_mod_cov.reshape( nbb,1)-
tab_xi0.reshape(nbb,1) ) ) )
In [6]:
t_min = np.unravel_index(tab_chi2_cov.argmin(), tab_chi2_cov.shape)
#Best OmBh^2, Best bias
print(r'Best Omegab*h^2: ',arrob[7])
print(r'Best bias:',tab_bias[50])
In [7]:
tab_delta_chi2 = tab_chi2_cov - tab_chi2_cov.min()
lev = [2.3,6.17,9.21]
plt.pcolormesh(arrob,tab_bias,np.transpose(tab_chi2_cov))
plt.contour(arrob,tab_bias,np.transpose(tab_delta_chi2),lev,linewidths=3,colors='white')
plt.plot([0.018896551724137931,0.018896551724137931],
[2.0101010101010104,2.0101010101010104],"o",color='red')
plt.xlabel(r'$\Omega_b h^2$', fontsize = 20)
plt.ylabel(r'Bias', fontsize = 20)
plt.show()